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Appendix  A 


1.0  INTRODUCTION 


This  report  is  the  final  report  on  Contract  No.  DAHC 
19-69-C-0032  with  the  Advanced  Research  Projects  Agency  entitled 
"Stress -’Wave  Propagation  Through  Earth-Water  Systems."  The  fun¬ 
damental  objective  of  this  study  is  to  develop  numerical  tech-, 
niques  to  treat  the  general  two-dimensional  stress  wave  propaga¬ 
tion  problem  through  nonlinear  earth  materials  including  the  ef¬ 
fects  of  water  flow  through  the  earth  materials. 

Prior  to  the  beginning  of  this  study,  a  numerical  tech¬ 
nique  was  developed  to  treat  the  dynamic  wave  problem  through 
arbitrary  nonlinear  media  (Ref.  5,  6  and  7)  without  including  the 
effects  of  water  on  the  propagation  process.  This  numerical  ap¬ 
proach  is  based  upon  the  finite  element  method  of  analysis  and 
led  to  the  development  of  a  large  computer  program  (termed  the ‘ 
SLAM  Code  for  identification,  the  acronym  standing  for  Stress 
Waves  in  Layered  Arbitrary  Media)  to  treat  either  the  general 
axisymmetric  or  plane  (stress  or  strain)  geometric  configuration. 
The  finite  element  approach  has  been  taken  in  this  development  to 
allow  the  user  a  general  flexibility  in  treating  two  dimensional 
problems  of  rather  complex  geometry  (inclusions,  material  layer¬ 
ing,  complex  boundaries,  etc.). 
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The  majority  of  computer  codes  that  have  been  developed 
to  treat  these  problems  are  based  upon  finite  difference  formula¬ 
tions  and  thus  would  have  difficulty  in  treating  complex  geometries 
of  interest.  In  the  finite  difference  formulations,  relatively 
uniform  meshes  are  distributed  over  the  half-space  of  interest  and 
equations  of  motions  applicable  at  each  node  are  developed  by  suit¬ 
ably  differencing  the  system  partial  differential  equations  of  mo¬ 
tion.  By  contrast,  in  the  finite  element  formulation,  the  two  di¬ 
mensional  half-space  is  divided  into  small,  arbitrarily  located, 
small  elements,  these  elements  being  connected  to  each  other  at 
their  vertices  or  riode  points.  From  the  element  type  and  geometry, 

equations  of  motion  are  developed  for  each  node.  The  primary  ad- 

•  % 

vantage  of  this  method  over  its  counterpart  difference  formulation 
thus  lies  in  its  flexibility  in  specifying  node  distributions  for 
a  particular  problem.  It  can  then  be  used  to  treat  problems  of  as 
complex  a  geometry  as  desired  with  no  change  in  computer  code 
formulation. 

Howe'er,  this  added  flexibility  is  obtained  at  a  cost, 
this  cost  being  the  increased  amount  of  information  that  must  be 
maintained  in  the  computer  (or  its  auxiliary  storage  units) .  Thus 
to  perform  a  calculation  for  a  simple  problem,  with,  say,  a  uni- 
form  mesh,  the  finite  element  formulation  will  be  inefficient  as 
compared  to  finite  differences.  For  problems  where  simple  mesh 
formations  can  be  used,  the  finite  element  solution  will  in  gen¬ 
eral  take  more  machine  time  to  achieve  the  same  results  as  the 
difference  solution.  Thus  the  very  flexibility  which  is  its 
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advantage  can  be  considered  a  disadvantage  for '-these  cases. 

After  the  development  of  SLAM  Code,  various  problems  of 
interest  were  investigated  to  determine  the  effects  of  material 
nonlinearities  on  the  wave  propagation  prooe:  s  (Ref.  5).  In  gen¬ 
eral,  three  types  of  problems  are  of  interest  wnen  studying  dy¬ 
namic  processes  through  earth  media,  in  the  first  type,  the  half 
space  is  subjected  to  high  intensity  pressure  loadings  caused  by 
high  energy  explosions.  The  resulting  ground  shock  effects  are 
highly  transient  and  are  characterized  by  relatively  short  dura¬ 
tion  shock  waves  of  high  strength.  Of  particular  interest  for 
this  proolem  is  the  rate  of  decay  of  the  shock  front  as  it  moves 
through  the  ground.  Clearly,  nonlinear  properties  of  the  mate¬ 
rial  significantly  influence  the  'decay  of  the  shock  strength 
since  large  non-recoverable  volume  changes  can  decrease  the  peak 
pressures  of  the  shock  front.  In  the  second  problem  type,  the 
half-space  is  subjected  to  ’ong  duration  low  intensity  vibratory 
type  loadings  associated  with  eartnquake  motions.  Again  nonlin¬ 
ear  properties  and  volume  changes  of  the  earth  material  signifi¬ 
cantly  alter  the  characteristics  of  the  motion  histories  sustained 
at  the  surface  of  the  ground. 

A  third  type,  not  often  considered  with  the  first  two 
problems,  is  concerned  with  the  response  of  near  shore  structures 
and  soils  subjected  to  ocean  storm  waves.  In  these  problems  non¬ 
linear  soil  properties  as  well  as  induced  seepage  and  pore  pres¬ 
sure  effects  tend  to  rduce  the  strengtn  of  the  soil  leading  to 
potential  massive  failures  of  the  soil  systems.  It  is  well  known 
(Ref.  11)  that  storm  waves  can  induce  slide  failures  on  even  flat 
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slopes  due  to  the  excess  seepage  conditions  developed  within  the 
ocean  sediments.  In  these  problems,  the  stress  and  motion  histo¬ 
ries  sustained  at  any  point  in  the  ground  are  significantly  influ¬ 
enced  by  the  nonlinear  characteristics  of  the  materials,  asso- 
ciated  volume  changes  and  any  induced  seepage  conditions  that  may 
occur.  Thus,  to  seriously  treat  the  time  dependent  response  of 
earth  media,  the  effects  of  pore  water  must  be  suitably  taken 
into  account. 

In  the  following  section  of  this  report,  a  brief  review 
of  the  analytic  approach  and  numerical  techniques  used  to  obtain  • 

i®  presented.  After  this  review,  the  results  obtained 
for  four  different  problem  types  are  presented,  these  problems 
being  chosen  primarily  to  indicate  the  areas  for  wnich  this  anal¬ 
ysis  is  applicable.  In- the  first  results  presented (Section  1.0), 
quasi-static  problems  are  considered;  that  is,  inertial  effects 
are  neglected.  The  problems  are  limited  to  the  (time-dependent) 
two  dimensional  consolidation  situations.  This  is  followed  (Sec¬ 
tion  4.0)  by  the  analysis  of  stress  wave  propagation  through  sat¬ 
urated  soil  columns  to  determine  the  influence  of  pore  fluid  on 
stress  wave  speed.  Sections  5.0  and  6.0  present  results  obtained 
for  the  effects  of  soil  liquefaction  generated  by  ocean  storm 
waves,  the  first  concerning  liquefaction  of  foundation  soils  and 
the  second  concerning  embankment  stability  to  the  storm  waves. 
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2 . 0  GOVERNING  SYSTEM  EQUATIONS 

In  Ref.  1,  the  derivation  of  the  system  of  equations 
governing  both  the  equilibrium  of  the  nodes  as  well  as  pore  water 
migration  (seepage)  was  presented  in  detail.  The  basic  assump¬ 
tions  made  are: 

(a)  the  soil  can  behave  as  a  nonhomogeneous,  aniso¬ 
tropic  and  inelastic  solid; 

(b)  the  pore  water  is  considered  as  incompressible 
with  respect  to  the  soil  component; 

(c)  the  deformation  of  the  soil  depends  upon  inter¬ 
granular  stresses  only; 

(d)  the  water  flows  through  the  soil  according  to 

.  I 

D'arcy’s  law. 

k 

For  completeness  of  this  report,  the  derivation  of  the  system 
equations  has  been  summarized  and  presented  in  Appendix  A.  The 
system  equations  can  be  written  symbolically  as 

W  * 

(i) 

[?}  ,  *•  LmI i*\ 

In  equations  1,  the  vectors  [ru]  and  {fw^  are  the  horizontal  (u- 
direction)  ana  vertical  (w-direction)  forces  generated  at  the 
nodes  of  the  finite  element  mesh.  The  matrices 

uu  uw  wu 

kt„,,  are  the  usual  elastic  stiffness  matrices  and  are  used  to  com- 

WW 9 

pute  the  elastic  components  of  the  forces  developed  at  the  node 
points  due  to  the  relative  displacements  of  the  element  nodes. 

The  vectors  ^u^  and  W  are  the  >v  r^zontal  and  vertical 
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displacements  of  the  nodes  of  the  mesh,  while  the  vector  {tt\  rep¬ 
resents  the  excess  pore  pressures  developed  at  the  node  points. 
The  matrices  £u  and  then  convert  those  excess  pore  pres¬ 

sure  into  equivalent  node  point  loadings. 

The  forces  at  che  node  points,  {fu}  and  Jfw^,  can  have 
three  components,  namely, 

(a)  applied  node  point  loads  due  to  applied  pressure 
loads ; 

(b)  fictitious  correction  forces  to  account  for  non¬ 
linear  stress-strain  behavior  of  the  soil  skeleton; 

(c)  inertia  loadings  for  dynamic  problems. 

The  forces  can  then  be  written  as 

W*  •  mm 

<•  {f2\  -  ImIM  <2) 

The  terms  with  the  superscript  A  in  equation  2  are  the  horizontal 
und  vertical  components  of  any  forces  applied  to  the  nodes  (from 
concentrated  loads  or  pressures  applied  to  specific  surfaces  in 
the  problem).  The  terms  with  the  superscript  N  refer  to  the  fic¬ 
titious  correction  forces  that  are  applied  to  cne  nodes  to  ac¬ 
count  for  any  nonlinearity  in  material  stress -strain  behavior  (or 
deviations  from  the  elastic  case) .  .  The  final  terms  of  equation  2 
refer  to  inertia  forces  if  dynamic  effects  a—  included. 

Al‘e  two  of  equation  1  then  represent  the  equilib¬ 

rium  of  total  strr-ses  at  a  point  in  the  half-space.  The  third 
equation  of  equat  un  1  represents  the  seepage  flow  relationship 
generated  from  D’arcy's  law.  The  matrix  [h}  is  dependent  upon 
e  coefficients  of  permeability  of  the  material  (as  well  as 
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properties  of  the  finite  element  configuration),  and  the  vectors 
iu\  and  ^w|  represent  the  nodal  velocities. 

The  vector  represents  the  rate  of  volume  change  of 
the  fluid  component  associated  with  each  node  point.  For  incom¬ 
pressible  fluid,  these  components  are  7,ero  for  interior  node 
points  (all  fluid  that  flows  into  an  element  must  flow  out) , 
while  for  some  boundary  node  points  (for  which  tho  excess  pore 
pressure  is  specified) ,  these  components  indicate  the  volume  of 
water  flowing  out  of  the  nodes. 

2 . 1  Solution  Procedures 

At  each  node  point  (except  at  boundary  nodes  where  ei¬ 
ther  displacements  and/or  excess -pore  pressures  are  specified), 
three  unknowns  must  be- determined  at  any  instant  of  time,  namel" 
the  two  node  displacements  (u  and  w)  and  the  excess  pore  pressure 
(tt)  •  The  solution  procedure  then  depends  upon  the  particular 
conditions  of  the  problem. 

For  example,  for  steady  state  flow  conditions  where  the 
coefficients  of  permeability  are  assumed  to  be  independent  of  tha 
intergranular  stress  (or  strain)  state,  the  time  variation  of  the 
parameters  is  zero.  Referring  to  equation  1,  this  allows  for  a 
direct  solution  of  the  third  equation  (by  elimination  procedures, 
say)  for  the  excess  pore  pressure  distribution.  With  these  as 
knowns,  the  first  two  of  equation  1  can  then  be  solved  for  the 
node  displacement  components.  If  the  material  is  elastic,  the 
nonlinear  correction  forces  of  equation  2  are  zero  so  that  a 
single  solution  is  all  that  is  needed.  If  the  material  stress- 
strain  relation  is  nonlinear,  the  displacement  solution  must  be 
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iterated  upon  since  the  nonlinear  correction  forces  are  functions 
of  the  node  displacements. 


For  problems  in  which  transient  effects  are  included 
but  for  which  inertia  effects  can  be  neglected  (slow  transients 
as  in  consolidation  problems) ,  the  following  procedure  can  be 
used.  The  velocity  intergration  procedure  is  based  upon  a  simple 
linear  velocity  approximation  during  a  small  time  step,  or 


V  v  /  *  9  v 


where  %i  represents  a  displacement  at  time  i,  *i_,  the  displace- 

•  • 

ment  at  the  preceding  time  3tep,  and  %L  represent  the  corre¬ 


sponding  velocities  and  At 


is  the  time  increment  between  t-i 


and  i.  Substituting  equation  3  into  equation  1  leads  to  the 
equations 


where  the  subscripts  i  refer  to  the  current  times.  The  "forcing" 
vector  can  be  written  as 

4  t^wV  l  ’"’t-i  ♦  *  u-im\  J 

and  is  defined  in  terms  of  the  previous  displacement  ar.d  velocity 
history.  The  "stiffness"  matrix  of  equation  4  is  symmetric  and 
usual  solution  procedures  can  be  used. 

Considering  an  elastic  porous  material,  at  a  particular 
instant  of  time,  the  effective  force  vector  of  equation  4  is 
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!  known,  since  the  applied  loads  are  specified  and  the  G-vectoi*  can 
be  computed  from  the  previous  time  step.  Obviously,  the  solution 
nust  start  from  some  time  when  the  initial  conditions  are  speci— 
,fied.  The  unknowns  (u,  w,tr)  at  the  cuiitm  *-ime  can  then  be  ob¬ 
tained  by,  say,  a  simple  elimination  technique.  The  solution  at 
the  following  time  step  can  then  be  obtained  using  the  current 
solution  as  input,  etc.  In  this  fashion,  the  solution  is  marched 
out  in  time. 

For  nonlinear  material  behavior,  this  procedure  has  to 
be  modified  since  the  nonlinear  correction  forces  of  equation  2 
are  also  functions  of  the  current  displacements.  To  overcome 
this  situation,  an  iteration  procedure  is  used  at  each  time  step 
to  ensure  that  the  nonlinear  correction  forces  are  properly  in¬ 
cluded  in  the  analysis;  At  each  time  step,  an  initial  trial  so¬ 
lution  is  obtained  by  using  approximate  values  for  these  correc¬ 
tion  forces  (usually  the  forces  from  the  previous  time  step) ,  and 
these  are  modified  by  the  iteration  process. 

For  problems  in  which  inertia  effects  are  no  longer 
negligible  (say,  stress  wave  problens) ,  a  similar  procedure  can 
be  used,  except  employing  a  linear  acceleration  extrapolation 
procedure ,  or 

**'  *l,  4  i  (*»-.  ♦  • 

*  *«-,  ♦  (6) 

Substituting  equation  G  into  equations  1  and  2  leads  to  another 
implicit  integrate  n  scheme  which  can  be  written  as 
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(7) 
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The  stitfness'  matrix  of  equation  7  i3  again  symmetric  and  may 
be  solved  by  s'mple  elimination  methods.  The  forcing  functions 
of  equations  7  are  defined  by 
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Again  the  solution  can  be  marched  out  in  time  as  before. 

A  draw  back  of  equation  8  occurs  for  cases  where  the 
time  increment  required  is  very  small  or  the  mass  vector  very 
large.  For  these  situations,  the  inertia  terms  of  the  "stiffness" 
matrix  may  become  extremely  large,  masking  out  the  stiffness 
terms.  This  in  turr.  may  lead  to  instability  in  the  implicit  in¬ 
tegration  procedure  3  used.  For  these  cases,  a  modified  approach 
can  be  used.  The  third  of  equation  1  can  be  solved  for  the  pore 
pressure  vector  (requiring  an  inversion  of  the  H  matrix)  and  this 
substituted  into  the  remaining  t.wo  force  equations.  For  this 
case,  the  force  equations  become 
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With  this  approach,  either  implicit  or  explicit  integration 
schemes  can  be  used. 
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3 . 0  SOIL  TEST  CONFIGURATIONS 

With  the  developed  computer  program,  numerical  results 
were  generated  for  several  soil  configurations  similar  to  the  usual 
soil  tests,  consolidation  and  triaxial  compression.  The  first  set 
of  data  assumed  elastic  soil  behavior,  since  for  these  problems 
analytic  solutions  are  available  or  can  be  easily  developed  for  com¬ 
parison  purposes.  The  first  nonlinear  soil  model  investigated  made 
use  of  a  Coulomb-Mohr  elastic  plastic  model  based  or.  the  concepts 
of  the  theory  of  plasticity. 

Although  this  model  is  often  used,  it  is  not  adequate  for 
modeling  stress-strain  behavior  (except  in  a  crude  sense)  and  would 
be  of  questionable  value  when  studying  pore  pressure  dependent  prob- 

*  t 

lems.  A  more  detailed  soil  model  was  then  investigated  which  ade¬ 
quately  predicts  stress  strain  behavior  of  a  particular  sand  sample 
and  was  developed  by  fitting  the  parameters  of  this  model  to  avail¬ 
able  experimental  data. 

3.1  One-Dimensional  Elastic  Consolidation 

The  first  problem  investigated  was,  naturally,  that  of 
the  classical  one-dimensional  consolidation  of  elastic  material. 

The  analytic  solution  available  for  comparison  is  the  standard 

e 

Terzaghi  solution  (Ref.  2).  The  problem  parameters  chosen  for  the 
investigation  are  shown  in  Figure  1.  The  computed  settlement-time 
history  at  the  top  of  the  soil  surface  is  shown  in  Fig.  2  and  com¬ 
parisons  made  with  the  exact  analytic  solution.  The  excess  pore 
pressures  developed  at  the  bottom  of  the  layer  are  shown  in  Figure  3 
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while  those  at  a  point  nearer  the  surface  are  shown  in  Figure  4. 
Again  comparisons  are  made  with  the  exact  solution  and  in  all  cases 
they  show  excellent  agreement.  The  pore  r ressure  distribution  at 
various  times  through  the  layer  is  shown  in  Figure  5.  Since  the 
pore  pressure  is  assumed  to  vary  linearly  within  a  given  element 
the  distribution  curves  are  piecewise  linear.  If  in  the  actual 
problem  the  pore  pressure  variation  is  share,  smaller  element  sizes 
must  be  used  to  suitably  approximate  the  solution. 

3 . 2  Triaxial  Elastic  Soil  Configuration 

The  second  model  considered  was  the  triaxial  soil  con¬ 
figuration  shown  in  Fig.  6a.  The  soil  model  was  considered  to  be 

t 

elastic  and  a  50  psi  vertical  pressure  applitd  at  the  initial  or 
zero  time.  The  finite  element  model  used  is  shown  in  Fig.  6b  and 
consists  of  28  rectangular  elements  to  represent  the  upper  quarter 
of  the  triaxial  sample.  The  elements  are  thus  axisymmetric  or  ring 
elements . 

To  obtain  the  analytic  solution,  it  was  assumed  that 
strain  conditions  in  the  sample  are  uniform.  The  initial  pore  pres¬ 
sure  developed  in  the  sample  (prior  to  drainage  occurring)  is  found 
from  the  following  analysis.  The  volume  change  per  unit  soil  volume 

9 

for  the  elastic  soil  is 

AV-  (10) 

where  the  barred  stresses  represent  the  intergranular  stresses,  and 
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E  and  V  are  the  elastic  modulus  and  Poisson's  ratio,  respect ive-y. 
Since  no  seepage  occurs  during  the  initial  conditions,  the  volume 
change  is  zero  or 


^  -  -  (Gy  -t  O'e') 


In  addition, 
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where  p  i9  the  excess  pore  pressure  and  CT^-  is  the  vertical  applied 
stress.  Combining  equations  11  and  12  leads  to  the  solution 

V4*/3  -  U3, 

The  initial  compression  of  the  soil  sample  is  simply 


T 


The  final  stresses  in  the  soil  system  are  obtained  when 
p  is  zero  (no  pore  pressure)  and 


*  ^AT 
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while  the  final  compression  of  the  soil  sample  Is 


A,  , 

4  E 


(16) 


where  L  is  half  the  original  sample  height*.  (height  of  the  finite 
element  model) .  The  settlement  from  the  initial  condition  to  the 
final  condition  is  governed  by  the  onc-dimensional  consolidation 
model  (since  one-dimensional  seepage  occurs  through  the  top  surface 
only,  with  the  modification  that  the  definition  of  the  coefficient 
of  consolidation  is 
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The  solution  to  the  particular  problem  of  Fig.  6  was  ob- 
tained  numerically  using  a  time  increment  of  0.1  seconds.  The  pore 
pressure  distribution  along  the  centerline  elements  is  shown  in 
Fig.  7  together  with  comparisons  with  the  analytic  solution.  As 
can  be  seen,  the  comparisons  are  excellent,  except  duriny  the  ear  Im¬ 
part  of  the  solution.  In  an  attempt  to  uncover  the  couse  of  the 
discrepancies,  the  same  problem  was  investigated  with  differing 

time  increments,  and  the  results  are  shown  in  Fig.  8.  As  may  be 

» 

noted  by  comparing  iige.  7  and  H ,  the  early  time  oscillations  found 
for  the  top  element  (Rlenumt  1'  are  related  to  the  time  step.  As 
the  rime  stop  is  decreased,  the  oscillations  disappear.  A  compar¬ 
ison  with  the  exact  solution  shows  that  the  computed  solution  is 
nliyhlly  lower  and  this  can  bo  nit.:  i  huted  to  the  fact  that  the  pore 


3.4 


12 


i 

ii 


i 


i| 

:j 


ii 

i| 


10 

*1 

ig 

rg 


io 

ft 


0 


pressure  profile  is  assumed  to  be  linear  across  the  element  while 

the  actual  pressure  profile  is  curaed,  particularly  at  the  early 

times . 

The  comparison  will,  the  middle  element  (Element  4)  is  not 
as  clear  cut,  however.  As  may  be  noted  from  Fig.  8,  this  pore  pres¬ 
sure  shows  an  initial  increase  in  pore  pressure  before  the  antici¬ 
pated  decay  occurs,  and  this  increase  is  independent  of  time  incre¬ 
ment  of  the  integration.  Since  this  phenomenon  did  not  occur  in 
the  elastic  plane  problem  discussed  previously,  it  must  be  concluded 
that  this  variation  is  concerned  with  the  coarseness  o £  the  finite 
element  mesh  in  the  radial  direction  for  this  axisymmetric  problem. 

NO  further  numerical  studies  have  been  conducted  on  this  problem 
as  yet,  however. 


3*3  Trlaxlal  CouJomb-Mohr  Model 

The  fust  triaxial  problem  including  nonlinear  material 
properties  that  was  investigated  was  the  same  model  shown  in  Fig.  6 
but  with  nonlinear  prop-rties  described  by  the  Coulomb-Mohr  yield 
condition  (Ref.  3).  For  stresses  within  the  yield  surface,  the 

soil  is  assumed  to  behave  elastically,  where  the  yield  surface  is 
defined  by 
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For  the  axisymmetric  stress  condition 


of  interest  for  this  problem, 
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where  the  bar  again  indicates  intergranular  stresses.  The  coeffi¬ 
cients  {oi,  k)  are  related  to  the  usual  strength  parameters  obtained 
from  a  tnaxial  tost  series,  4>  ,  the  angle  of  internal  friction,  and 
c,  the  cohesion,  by 
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For  stresses  on  the  yield  surface,  plastic  strain  components  are 
determined  fron  the  usual  normality  principal. 

Prior  to  investigating  this  problem  numerically,  the  ana¬ 
lytic  solution  for  the  initial  stress  condition  was  obtained  (no 
drainage  allowed).  As  the  vertical  stress  is  slowly  increased,  the 
soil  behaves  elastically  and  the  previous  solution  applies.  Sub- 
at.  tuting  equations  12  and  13  into  equations  19  yields 
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For  plastic  yielding  to  begm,  the  critical  vertical  stress  must 
roach  the  vnluo 
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tete. 


Fo’-  applied  stresses  larger  than  this  critical  value,  plastic  fl  w 
must  be  accounted  for,  making  use  of  the  normality  relation 
(Ref,  3) ,  which  for  this  problem  becomes 
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where  (  £r  ,  )  are  the  radial  and  vertical  components  of  the  pla 

tic  strain  rate  vector.  The  plastic  volume  change  iB 

•  i 
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where  is  the  total  plastic  vertical  strain,  while  the  elastic 
'Glume  change  is 


Knowing  that  the  total  volume  change  is  zero  (no  drainage  out  of 

t 

the  sample  is  allowed) ,  the  solution  can  be  readily  obtained  for 
any  applied  stresses  greater  than  the  critical,  or 


The  results  for  a  particular  undrained  case  are  shown  in 
Fig.  9.  The  vertical  pressure  is  applied  "slowly"  with  a  rise  time 
of  50  seconds  until  it  reaches  a  peak  pressure  of  50  psi.  The  par¬ 
ticular  properties  of  the  soil  chosen  were 

F  -  1  0  CO  ys  iL 

l)  -  0.25 
e  -  2o,  8  psi 

* 

<$)  -  3 0° 

For  this  condition  the  critical  vertical  stress  is  reached  when  '3'1_ 

rsj 

is  43.4  psi  and  the  corresponding  pore  pressure  is  14.45  psi.  As 
the  vertical  stress  is  increased  to  50  psi,  plastic  flow  takes  place 
(along  with  plastic  volume  exDansion)  and  the  pore  pressure  reduces 
to  11.1  psi.  Five  computer  runs  were  made  for  this  problem  using 
different  time  Bteps  as  seen  in  Fkj.  9.  In  each  case,  the  nonlinear 


3.8 


i 


correction  forces  in  the  equilibrium  equations  ..were  taken  as  the 

»• 

|  value  computed  during  the  previous  time  step.  As  can  be  noted, 

the  smaller  the  time  step,  the  better  the  approximation,  as  expected. 
As  an  alternate  to  this  procedure,  the  nonlinear  correction  forces 
in  a  given  time  step  can  be  recomputed  by  iteration  (obtain  a  trial 
solution,  computed  correction  force,  obtain  new  solution,  etc.). 

For  this  problem  of  proportional  loading,  this  procedure  is  equiva¬ 
lent  to  using  smaller  time  steps  without  iteration  during  each  time 
step. 

After  the  final  equilibrium  condition  is  reached  under  no 
drainage  conditions,  the  drained  sit  uation  can  be  achieved  by  let¬ 
ting  the  pore  pressure  decrease  to  zero  by  allowing  drainage  through 
the  top  and  bottom  surfaces  of  the  soil  sample.  It  can  be  shown 
that  for  this  soil  model,  the  decay  of  the  pore  pressure  will  occur 
•lA*tically;  that  is,  the  intergranular  stress  state  will  move  off 
the  yield  surface  as  the  pore  pressure  decreases,  so  that  the  decay 
rate  will  be  as  described  in  the  previous  elastic  triaxial  solution. 

The  solutions  for  these  cases  are  shown  in  Fig.  10  where 
the  vertical  intergranular  stress  is  plotted  as  a  function  of  the 
total  vertical  strain  for  various  values  of  the  cohesion  and  for  a 
fixed  value  of  the  friction  angle  of  30*.  If  t  .  cohesion  is  24.0 
psi  or  greater,  the  soil  sample  always  remains  elastic.  The  ini¬ 
tial  stress  state  when  a  vertical  stress  of  50  psi  is  applied  and 
no  drainage  is  allowed  is  CT^  »  33.3  psi  and  p  »  16.7  psi.  When 
drainage  is  then  allowed,  the  pore  pressure  decreases  to  zero,  and 
the  vertical  intergranular  stress  increases  to  the  applied  stress 
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of  50  pal.  The  final  state  is  tho  same  a:  would  occur  if  the  sample 
had  been  tested  dry  (no  pore  pressure).  If,  however,  the  cohesion 
is  lower ,  the  initial  undrained  state  causes  plastic  flow  to  ooour 
with  the  drained  condition  occurring  elastically,  as  shown.  For  a 


value  of  cohesion  equal  to  14.4  psi,  the  initial  undrained  state 
occurs  with  no  excess  pore  pressure,  and  the  drained  state  is  the 
same  as  the  undrained  state.  For  values  of  cohesion  less  than  14.4 
psi,  equilibrium  under  the  applied  loads  cannot  be  maintained.  It 
should  be  pointed  out  that  for  values  of  cohesion  between  14.4  and 


24,0  pni,  the  dry  test  will  show  no  plastic  flow,  while  the  undrained** 

•  i 

drained  sequence  will  yield  plastic  strains. 

It  is  clear  then  that  even  for  this  relatively  simple 
soil  model,  the  stress-strain  behavior  between  saturated  and  unsat¬ 
urated  soil  samples  will  be  different  and  will  be  influenced  by 
the  rate  of  loading  (as  compared  with  the  rate  of  pore  pressure 
decay).  To  investigate  this  analytic  solution  further,  the  pre¬ 
vious  solution  was  nondimensionalized  in  the  following  fashion. 

Non  dimensional  parameters  are  defined  as 
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The  upper  and  lower  limits  of  cohe.'sion  for  which  a  nonlinonr  solu¬ 
tion  (stable  plastic  strains  will  occur)  can  be  obtained  for  t.he 

undrained  case  are 
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For  any  value  of  cohesion  between  these  limits,  the  solution 


yields 
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After  this  initial  solution  occurs,  the  addition  vertical  strain 
that  will  develop  at  the  excess  pore  pressure  is  allowed  to  decay 


to  zero  is 


so  that  the  final  strain  is  the  sum  of  the  strains  from  equations 


3. 


28  and  29.  The  solutions  for  several  parameter  variations  are  shown 

in  Figs.  11  through  14.  In  Fig.  11,  the  nondimensional  vertical 

intergranular  stress  is  plotted  as  a  function  of  the  ratio 

D  *  * 

where  is  the  vertical  strain  that  would  occur  in  the  dry  state 
and  is  simply 


(30) 


As  may  be  noted,  the  difference  in  limiting  values  of  cohesion  for 
this  problem  is  relatively  small,  but  the  influence  on  the  final 

i 

strain  is  large  (ratio  of  6.25).  Curves  are  shown  for  four  equally 
spaced  values  of  cohesion  between  the  limiting  values. 

The  same  solution  is  shown  in  Fig.  12,  except  that  the 
friction  angle  was  increased  from  5°  to  30°.  As  may  be  noted,  the 
final  strains  are  much  lower  than  those  of  Fig.  11,  and  the  associ¬ 
ated  plastic  strains  occurring  during  the  initial  undrained  state 
are  much  smaller.  This  is  due  to  the  fact  that  for  the  higher  fric¬ 
tion  angle  the  plastic  volume  expansion  is  larger  than  for  the 
smaller  friction  angle  causing  the  excess  pore  pressure  to  decay 

9 

more  rapidly  as  plastic  strains  develop.  Fig.  13  shows  the  same 
results  for  a  still  larger  friction  angle  of  45°,  again  showing  a 
smaller  difference  in  final  strains. 

The  results  for  a  different  value  of  Poisson's  ratio 
(Z^=  0.25)  are  shown  in  Fig.  14  fur  a  friction  angle  of  30°.  Ah 
can  be  noted,  the  behavior  is  essentially  different  than  that  of 
Fig.  12.  This  is  due  to  the  fact  that  the  elastic  volume  change 
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during  the  initial  loading  decreases  as  Poisson's  ratio'  increases. 

3.4  McCormick  Ranch  Sand  Model 

It  is,  of  course,  well  known  that  the  relatively  simpli¬ 
fied  constitutive  models,  such  as  the  Coulomb-Mohr  model,  can  only 
crudely  approximate  the  stress-strain  behavior  of  real  soils.  In 
order  to  properly  take  into  account  the  influence  of  pore  fluid  on 
sell  response,  more  realistic  models  must  be  developed.  An  example 
of  such  a  model  was  presented  in  Ref.  4  wherein  the  parameters  of 
the  model  were  chosen  to  match  (as  closely  as  possible)  available 

•  4 

experimental  data  on  a  particular  sand  sample,  known  as  McCormick 
Ranch  Sand.  A  rather  extensive  series  of  triaxial,  uniaxial,  and 
hydrostatic  compression  tests  were  conducted  and  an  attempt  was 
made  to  fit  the  analytic  model  so  as  to  reproduce  the  available  data. 

It  was  found  that  for  the  particular  parameters  chosen 
the  stress-strain  curve  during  the  initial  load-unload  cycle  could 
be  adequately  reproduced  for  the  triaxial  compression  test  (over  a 
wide  range  of  confining  pressures)  and  for  the  uniaxial  compression 
test.  The  soil  model,  however,  was  significantly  stiffer  under 
hydrostatic  compression  (although  the  shape  of  the  load-unload  curve 
was  the  same)  than  the  experimental  data. 

The  model  is  based  on  the  following  analysis.  The  hydro¬ 
static  and  deviatorio  stress-strain  components  are  related  by 
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where  ^  1b  a  critical  hydrostatic  pressure  (positive  in  compres¬ 
sion  and  Jj[  is  the  second  invariant  of  the  deviatoric  stresses 
(equation  19) . 

To  match  the  specific  test  results  for  the  sand  sample, 
the  following  parameters  were  found  to  best  reproduce  all  the  data: 
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The  stress-strain  behavior  under  uniaxial  compression  is 
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shown  in  Fig.  15  under  both  initial  loading  conditions  as  well  as 
strain  load/unload  cycling.  As  may  be  noted,  the  stress-strain 
response  exhibits  the  characteristic  stiffening  effect  as  well  as 
the  nonrecoverable  behavior  under  load  cycling.  The  pressure  ranges 
shown  are  higher  than  normally  used  but  suitable  modification  of 
the  data  input  would  convert  this  typical  response  to  lower  stress 
ranges  of  interest. 

The  behavior  under  triaxial  compression  is  presented  in 
Figs.  16  to  19  and  again  exhibits  much  of  the  characteristics  antic¬ 
ipated  for  a  sand  sample.  During  the  load/unload  cycling,  the 
model  can  be  further  improved  to  reproduce  test  data  by  modifying 
the  shear  modulus  formulation  under  reload  conditions  to  better 
match  strain  behavior  with  constant  load  cycling. 

The  previous  data  were  obtained  for  the  Ranch  Sand  model 
in  the  dry  condition.  To  determine  the  behavior  with  pore  fluid, 
similar  problems  were  investigated  including  load  cycling  effects. 

In  Fig.  20,  the  triaxial  response  is  presented  for  a  consolidated/ 
undrained  experiment  with  load  cycling  m  the  vertical  direction  x 
corresponding  to  the  load  cycles  shown  in  Fig.  18  for  the  dry  sam¬ 
ple.  In  both  cases,  lateral  or  confining  stresses  were  maintained 
constant.  As  can  be  seen  in  Fig.  20,  the  effect  of  pore  pressure 
is  to  decrease  the  axial  strain  increment  between  load  cycles. 

That  is,  in  the  undrainod  state,  the  soil  model  "shakes  clown"  to 
effectively  a  linear  model,  although  strong  nonlinear  behavior 

again  takes  effect  as  the  applied  load  is  finally  increased  beyond 
the  load  cycling  regime. 
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Similar  behavior  is  shown  .in  Fig.  21  &here  the  applied 
vertical  load  is  cycled  through  tho  complete  load  range  from  0  to 
300  psi.  This  tesc  corresponds  to  the  dry  triaxial  test  shown  in 
Fig.  19.  Again,  it  may  be  noted  that  within  a  load  cycle,  pore 
pressure  effects  cause  the  stress-r.train  behavior  to  "shake  down" 
to  an  effective  elastic  state.  Of  course  this  type  of  response 
can  be  modified  by  changing  the  definition  of  the  reload  shear  mod¬ 
ulus  as  defined  by  equation  34.  A  plot  of  the  invariants  of  effec¬ 
tive  stresses  during  the  loading  cycle  for  the  triaxial  tests  is 
shown  in  Fig.  22,  for  both  the  consolidated  undrained  and  drained 
tests.  As  may  be  noted,  is  constant  during  the  undrained  test 
indicating  that  the  bulk  modulus  (equation  33)  is  constant  with  this 

•  s 

model  Therefore  the  cycling  response  will  be  completely  dependent 
upon  the  variation  in  the  deviatoric  response,  or  the  shear  modulus 
behavior.  The  cycling  response  will  be  eBBentially  elastic  as  long 
as  the  shear  modulus  is  maintained  as  the  unloading  modulus  within 
the  cycling  load  range. 

Two  other  triaxial  experiments  were  conducted  where  the 
samples  were  consolidated  under  a  confining  stress  of  400  psi, 
loaded  vertically  in  the  drained  state  to  630  psi  and  then  further 
loaded  cycled  between  575  psi  and  690  psi  in  both  the  drained  and 
undrained  states.  A  comparison  of  the  results  is  shoton  in  Fig.  23, 
in  which  anticipated  responses  were  determined. 
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Fig.  9  Triaxial  Compression,  Cculomb-Mohr  Material,  Undrained  Test 
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Fig.  14  Triaxial  Corapresoion  Test,  Coulonb-Xohr  Material,  Non-Dimensional  Solution 
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II 
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^  •  0  Stress  Wave  Propagation  in  Uniaxial  Compression 

The  first  dynamic  problem  considered  was  the  one¬ 
dimensional  soil  column  subjected  to  a  step  pulse  pressure  at 
its  end,  a  configuration  often  used  in  determining  wave  pro¬ 
pagation  properties  through  actual  soil  columns.  The  configu¬ 
ration  studied  assumed  the  soil  column  to  be  contained  within 
an  imprevious  membrane  so  that  no  pore  fluid  (for  the  case  of 
saturated  soils)  could  escape.  He  ever,  it  was  assumed  that 
the  membrane  provided  no  lateral  restraint  to  the  soil  column,, 
so  that  even  for  an  incompressible  fluid,  a  finite  compressed 
wave  speed  exists.  If  lateral  inertia  effects  are  neglected, 

the  effective  modulus  of  the  soil  (plane  strain/stress  condition) 
is 


E*  •  E  /  ( i  -  it l) 


(36) 


and  the  effective  dry  wave  speed  is 


.  ,  ,Vt  ' 


(37) 


The  finite  element  configuration  used  for  this 
problem  is  shown  in  Fig.  24,  and  consisted  of  a  set  of  square 
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plane  strain  elements  with  no  lateral  restraint  provided  to 
the  nodes.  In  the  numerical  calculation,  lateral  inertia 
effects  are  in  fact  included  since  each  node  is  considered 
to  have  two  degrees  of  freedom.  The  first  problem  run  con¬ 
tained  no  pore  fluid,  and  the  computed  stress  profile  along 
the  rod  is  shown  in  Figure  25.  As  can  be  noted,  the  computed 
stresses  show  the  usual  oscillations  about  the  true  solution 
caused  by  the  sharp  stress  front  applied  to  the  rod.  The 
propagation  of  the  pulse  through  the  soil  is  shown  in  Fig.  26 
and  the  computed  wave  speed  is  seen  to  compare  favorably  with 
the  actual  effective  wave  speed.  A  similar  result  is  shown  in 
Fig.  27  for  the  case  of  a  different  Poisson's  ratio. 

For  the  case  of  saturated  soils,  the  stress  profile 

i 

along  the  rod  at  a  late  time  in  the  pulse  is  shown  in  Fig.  28. 
As  may  be  noted,  the  wave  speed  is  higher  than  the  dry  wave 
speed,  since  the  pulse  is  already  reflecting  off  the  bottom 
boundary.  At  the  loaded  end  of  the  soil  column,  the  solution 
approaches  the  static  solution,  which  for  an  incompressible 
pore  fluid,  indicates  that  the  pore  pressure  is 


^  *  CT*  /  2. 


(38) 
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where  0^  is  the  applied  static  loading.  Figure  29  shows 
similar  results  for  undrained  soil  with  a  different  Poisson's 
ratio.  As  may  be  noted  by  comparing  Fxys.  28  and  29,  the 
effective  wave  speed  decreases  as  Poisson's  ratio  increases, 
which  is  opposite  to  the  results  computed  for  drv  soils. 

Obviously  the  reason  for  this  discrepancy  lies  in  the  amount 
of  lateral  inertia  effects  that  are  allowed  to  occur  for  the 
different  Poisson's  ratios. 

A  plot  of  peak  stress  traverse  through  the  soil 
rod  is  shown  in  Figure  30,  for  cases  of  both  dry  and  saturated 
soils.  As  may  be  seen,  the  effective  dry  wave  velocity  increases 
as  Poisson's  ratio  increases  (as  predicted  by  equation  37). 
However,  the  opposite  effect  occurs  for  saturated  soils  since  • 
the  partition  of  energy  into  t\\e  lateral  direction  essentially 
lowers  the  effective  wave  speed. 

A  series  of  problems  were  investigated  for  different 
coefficients  of  permeability.  Strep*:  profiles  along  the  rod 
are  shown  in  Fig.  31  for  various  valuer  of  K  where 


\  «  k/fcur 


(39) 


and  k  is  the  coefficient  of  permeability.  From  these 

# 

results,  it  may  be  determined  that  the  effective  wave  speed 
(as  determined  by  peak  stress  traverse  through  the  rod)  is 
unaffected  by  the  coefficient  of  permeability.  However,  the 
stress  profile,  both  in  the  vicinity  of  the  front  and  behind 
the  front,  is  affected.  This  would  indicate  that  local  fluid 
migration,  controlled  by  the  coefficient  of  permeability, 
directly  inf luence* the  local  intergranular  stress  conditions. 
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5 . 0  POTENTIAL  FOR  WATER  WAVE  INDUCED  SOIL  LIQUEFACTION 


A  problem  of  considerable  interest  for  off-shore  struc¬ 
tures  is  the  potential  for  liquefaction  of  soils  due  to  storm 
waves  traversing  a  site.  If  foundation  soils  liquefy,  then  ob¬ 
viously  they  lose  strength  and  ability  to  support  any  applied 
loads.  It  is  well  known  (Ref.  11)  that  pressure  changes  on  the 
sea  floor,  associated  with  passage  of  an  ocean  wave,  can  cause 
shear  fiilure  in  soft  sediments,  this  shear  failure  developing 
from  the  effects  of  both  water  flow  and  reduced  strength.  In  ac¬ 
tual  cases,  the  waves  impose  an  oscillatory  motion  on  the  soft 
sediments,  leading'  to  repeated  reversals  of  shear  strength,  re¬ 
molding  Ox  the  soil  with  an  associated  loss  of  strength. 

For  this  study,  an  approach  has  been  taken,  based  upon 
the  combined  stress-seepage  analysis  presented  in  Section  2.0, 
which  may  be  used  to  establish  the  potential  for  liquefaction  be¬ 
ing  induced  in  marine  soils  due  to  ocean  storm  waves.  The  con¬ 
figuration  analyzed  is  shown  in  Figure  32,  and  consists  of  a  half¬ 
space  of  submerged  saturated  soils  to  which  a  differential  total 
stress  is  applied  at  the  surface  equal  to  the  pressure  differen¬ 
tials  generated  by  the  storm  wave.  In  addition,  an  excess  pore 

pressure  distribution  along  the  surface  imposes  a  seepage 

# 

condition  aggravating  the  intergranular  stress  picture  in  the 
soil  due  to  the  seepage  forces  generated.  In  fact,  <_nis  storm 
wave  is  traversing  the  site  at  a  (say,  constant)  velocity  so  that 
the  storm  wave  pressure  distributions  at  any  one  location  are 
functions  of  time. 
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An  approximate  approach  to  the  problem  is  to  freeze 
the  wave  at  any  instant  of  time  and  analyze  the  foundation 
response  to  steady  load  and  seepage  criteria.  This  approxi¬ 
mation  neglects  the  potential  time  lag  effects  which  may  oc¬ 
cur  in  the  foundation  soil,  but  this  effect  cannot  be  too 
significant  for  pervious  soils  at  the  surface  of  the  founda¬ 
tion.  If  this  approach  is  utilized  and  if  the  soils  are  as¬ 
sumed  to  behave  linearly,  the  infinite  half-space  problem  can 
be  reduced  by  required  symmetry  conditions  to  analyzing  the 
restricted  region  shown  in  Figure  32. 

Applied  Pressure  Condition 

i 

As  mentioned  above,  the  loading  condition  applied 
to  the  soil  surface  consists  of  equal  amounts  of  total  pres¬ 
sures  and  excess  pore  pressures.  The  excess  pore  pressures 
will  induce  seepage  to  occur  from  the  area  of  the  wave  peaX 
to  the  wave  trough  and  the  equal  applied  surface  pressures 
are  applied  to  maintain  proper  intergranular  stress  condi¬ 
tions.  This  last  statement  may  be  amplified  by  considering  a 
uniform  increase  in  the  water  level  across  the  site,  leading 
to  a  uniform  increase  in  pore  pressure  along  the  soil  surface. 
No  seepage  occurs  for  ;his  condition,  but  the  intergranular 
stress  condition  will  change  unless  an  equal  surface  loading 
is  also  applied. 

As  a  wave  passes  a  site,  an  increase  in  pressure,  p, 
occurs  below  the  cres.,  while  beneath  the  trough  there  is  a 
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pressure  decrease,  -p.  The  characteristics  of  the  applied 
pressure  pulse  used  in  the  analysis  was  taken  as  a  sine  func¬ 
tion  of  the  form 


•y  -  p0  Stvt  ( 2»tx/d} 


where  p  is  the  pressure  increase  at  a  given  point,  Po  is  the 
peak  water  pressure,  and  D  is  the  wave  length  or  distance  be¬ 
tween  wave  peaks.  At  any  point,  the  corresponding  presbure 
var  at ion  with  time  is 


-p  *  P.  Cv*  *  /t) 


where  T  is  the  period  of  the  wave.  The  relation  between  D 


and  T  is  specified  by 


D  *  V  T 


where  V  is  the  velocity  of  the  wave  profile.  The  magnitude 
of  the  peak  pressure  change,  PQ,  which  is  in  phase  with  the 
wave,  depends  on  the  wave  length,  the  depth  of  water  and  the 
wave  height.  The  actual  pressure  puls**  "''olied  by  the  ocean 
wave  is  complex,  but  using  tables  ba>:  .ri  uj>cv  linear  wave 
theory  (Ref.  12)  the  following  parameters  vre  used 

P.  •  400  ptf 
V  *  V?  fps 

T*  m  (43) 

X)  «  sic  ft 
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These  correspond  to  a  storm  wave  of  about  30  feet  occurring  in  a 
water  depth  of  about  50  feet.  These  parameters  correspond  to  de¬ 
sign  conditions  for  near  shore  structures  which  may  be  placed  on 
the  East  Coast. 

5.2  Problems  Investigated 

Returning  to  Figure  32,  the  finite  element  mesh  is  used 
only  in  the  region  extending  across  one  side  of  the  loading  wave, 
since  by  symmetry  conditions  the  response  on  the  other  side  of 
the  loading  wave  can  be  deduced.  For  each  soil  site  considered, 
two  computer  runs  were  generated,  the  first  to  obtain  the  in-situ 
stress  conditions  due  to  the  bouyant  unit  weight  of  the  soil  and 
the  second  to  determine  the  additional  stresses  due  to  the  com¬ 
bined  effects  of  the  excess  pore  pressures  (seepage  stress  ef¬ 
fects)  and  applied  surface  pressures.  The  final  output  is  then 
the  total  stress  state  due  to  both  the  initial  and  final  load 
conditions.  For  inelastic  soils,  the  load  cycling  effect  could 
be  estimated  by  alternating  the  sign  of  the  applied  wave  loadings 
to  determine  growth  of  soil  strain  with  time,  while  for  elastic 
materials,  this  is  not  required. 

For  each  phase  of  the  problem,  the  specified  boundary 
conditions  used  for  the  finite  element  mesh  are  as  shown  in  Fig¬ 
ure  33.  For  each  problem,  these  conditions  satisfy  the  required 
symmetry  conditions  imposed  by  freezing  the  wave  motion.  Three 
separate  soil  site  configurations  were  investigated  for  this  prob¬ 
lem,  in  an  attempt  to  estimate  general  liquefaction  behavior,  and 
these  sites  are  shown  in  Figure  34.  The  first  site  consists  of  a 
pervious  sand  stratum  overlying  a  relatively  impervious  clay 
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stratum,  the  second  site  consists  of  a  sand/clay/sand  configura¬ 
tion,  and  the  third  is  a  clay/sand/clay/sand  site,  these  being 
typical  of  offshore  sites.  The  loose  surface  sands  are  relatively 
mobile  in  that  scouring  offer!  transport  these  soils  relatively 
easily. 

For  the  problems  investigated,  the  coefficients  of  per¬ 
meability  of  the  sands  was  taken  to  be  about  200  times  the  coef¬ 
ficients  of  the  clay  strata»  In  addition,  for  each  soil  layer, 
the  horizontal  coefficient  of  permeability  was  assumed  to  be  5 
tirr.es  larger  than  the  vertical  coefficient,  these  values  again 
beina  rather  typical  for  these  soils.  For  sites  1  and  2,  with 
the  surface  sand  layer,  it  can  be  anticipated  that  the  majority 
of  the  flow  pattern  will  take  place  through  the  sands,  leading  to 
large  seepage  pressures  through  this  zone.  In  addition,  for  site 
2,  a  variation  of  parameter  study  was  made,  varying  the  perme¬ 
ability  ratios  for  the  sand/clay  strata  to  determine  their  effect 
on  liquefaction  potential. 

5 . 3  Definition  of  Liquefaction  Fotential 

For  any  problem  investigated,  states  of  intergranular 
stress  and  excess  pore  pressure  distribution  can  be  determined. 

The  problem  remaining  is  to  interpret  these  stress  data  to  deter¬ 
mine  if  any  portion  of  the  soil  is  at  or  near  its  peak  shear 
strength.  If  the  computed  intergranular  stress  states  are  rela¬ 
tively  small  in  comparison  to  the  available  soil  strength,  it  can 
be  anticipated  that  excess  strength  is  available  and  no  liquefac¬ 
tion  could  take  place.  If  the  stress  states  exceed  the  available 
strength,  liquefaction  obviously  occurs  and  no  additional 
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strength  is  available  in  general  to  support  additional  structural 
loads . 


A  typical  way  of  defining  available  strength  is  to  use 
the  Coulomb-Mohr  failure  criteria  defined  by  equation  18,  or 

a  J,  *  /jy  --  k. 

(44) 

where  J,  is  the  hydrostatic  stress  invariant,  j'  is  the  devia- 
toric  stress  invariant  (effective  shear  strength)  and  the  parame¬ 
ters  ot  and  k  are  related  to  the  usual  strength  parameters,  C  and 
0,  obtained  from  a  triaxial  test  series,  and  are  defined  by  equa¬ 
tion  20.  This  form  of  the  strength  law  leads  to  the  usual  linear 
Mohr  strength  envelope  used  in  soil  testing  and  analysis.  If  an 
intergranular  stress  state  is  specified  at  a  point,  the  hydro¬ 
static  invariant  can  be  computed,  and  the  "available  shear 
strength"  defined  by 


'  ft'  ) AVAIL  r  fe.  -  *  J| 


(45) 


Comparing  this  value  with  the  actual  applied  deviator  stress, 
(/jJ"  )AppLiED*  one  can  assess  if  the  stress  state  is  at' or  near 
failure.  An  obvious  factor  of  safety  can  then  be  defined  as 


‘•Pi^'Ll  /Cfiv); 


applied 


(46) 


If  this  factor  is  near  unity,  it  can  be  anticipated  that  no  sig¬ 
nificant  strength  is  available  to  support  structural  loads,  and 
liquefaction  has  occurred. 


5.4  Numerical  Results 


The  finite  element  mesh  used  for  Site  No.  1  is  shown  in 
Figure  35,  and  consists  of  a  relatively  uniform  mesh.  The  excess 
pore  pressure  distribution  pattern  is  shown  in  Figure  36,  and  in¬ 
dicates  the  magnitude  of  the  pore  pressures  developed  by  the 
ocean  wave  above  the  ambient  pore  pressure  condition.  Obviously, 
by  symmetry,  the  pore  pressures  pattern  on  the  other  side  of  the 
right  boundary  of  the  mesh  is  just  the  mirror  image  of  the  pat¬ 
tern  shown.  The  horizontal  displacement  profile  under  the  posi¬ 
tive  wave  is  shown  in  Figure  37  and  generally  indicates  that  un¬ 
der  this  portion,  the  horizontal  stress  increments  developed  by 
the  wave  are  tensile  and  tend  to  reduce  the  confinement  of  the 
soil.  Under  the  negative  portion  of  the  wave,  the  opposite  would 
be  true. 

The  corresponding  vertical  stress  increment  due  to  the 
storm  wave  is  shown  in  Figure  38  and  as  may  be  noted  are  rela¬ 
tively  small.  Similarly  shear  stresses  (Fig.  39)  are  small.  The 
major  effect  is  then  to  generate  lateral  tensile  stresres  under 
the  crest  of  the  wave,  as  shown  in  Figure  42,  in  the  top  sand 
layer.  A  comparison  of  available  strength  with  applied  effective 
shear  stress  is  shown  in  Figure  43,  and,  as  may  be  noted,  near  the 
surface  of  the  sand  layer,  the  available  strength  approaches  the 
used  strength  indicating  a  low  factor  of  safety. 

Similar  data  for  Site  No.  2  is  shown  in  Figures  44  to 
48,  with  the  same  general  properties  of  the  stress  state.  Near 
the  surface  of  the  top  sand  layer,  the  safety  factor  approaches 
1.0,  and  in  fact  under  the  crest  of  the  wave  (left  boundary)  the 
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safety  factor  becomes  less  than  1.0.  For  Site  No.  3  (Figs.  49 
to  54)  ,  similar  resalts  are  obtained  with  the  exception  that  the 
safety  factors  in  the  top  clay  and  sand  layers  are  higher  than 
those  in  the  previous  problem.:.  This  indicates  that  the  pore 
pressure  dissipation  through  the  top  clay  layer  is  sufficient  to 
reduce  the  pore  pressures  and  thus  radial  stress  increments  in 
the  sand  layer  which  is  most  susceptible  to  the  pore  pressure 
effects. 

As  mentioned  previously,  a  variation  of  parameter  study 
was  made  for  Site  No.  2,  varying  the  both  the  ratio  of  sand  to 
clay  permeabilities  as  well  as  the  ratio  of  horizontal  to  verti¬ 
cal  permeabilities.  In  Figure  55,  the  safety  factor  profile  un¬ 
der  the  crest  of  the  wave  (left  boundary)  is  shown  as  a  function 
of  depth  as  the  ratio  of  sand  to  clay  permeability  was  varied 
from  1  to  200.  In  the  cases  investigated,  this  crest  boundary 
led  to  the  lowest  factor  of  safety.  As  may  be  noted,  the  safety 
factors  decreased  as  the  ratio  increased  (more  flow  restricted  to 
the  top  sand  layer)  .  For  this  problem,  the  entire  surface  sand 
layer  liquefies  (or  loses  all  its  strength) .  Figure  56  shows 
similar  data  for  a  variation  of  horizontal  to  vertical  permeabil¬ 
ity  in  each  layer.  As  may  b2  noted,  this  parameter  does  not  sig¬ 
nificantly  alter  the  results.  ♦ 
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Fig.  34-  Site  Configurations  Investigated 
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6 . 0  EMBANKMENT  STABILITY  TO  OCEAN  WAVES 

In  t  .is  section  of  the  report,  the  stability  analysis 
of  an  off-shore  embankment  is  presented,  this  embankment  being 
subjected  to  an  ocean  storm  wave  similar  to  that  used  in  the  pre¬ 
vious  problem.  The  profile  of  the  embankment  is  shown  in  Figure 
57  and  consists  of  a  concrete  caisson  protected  on  the  seaward 
side  by  an  embankment  of  coarse  material  varying  from  gravel  to 
large  "dolose"  pieces  on  the  external  side.  This  embankment  is 
placed  atop  a  foundation  similar  to  Site  No.  2  of  Section  5.0  and 
subjected  to  a  large  storm  wave  approaching  from  sea  (right  of 
Figure  57) .  The  objective,  of  course,  is  to  determine  the  ade¬ 
quacy  of  the  embankment  in  protecting  the  caisson  and  landward 
zone  within  the  embankment,  by  estimating  the  safety  factor 
against  a  potential  outward  slide  caused  by  the  wave. 

From  Section  5.0,  it  was  determined  that  the  surface 

sand  layer  will  liquefy  under  the  crest  of  the  wave,  so  that  as 

the  crest  approaches  the  toe  of  the  embankment,  a  potential  slide 

may  occur  since  the  sand  layer  has  little  available  shear  strength 

to  support  the  embankment.  Since  it  was  anticipated  that  a  po-  * 

tential  failure  condition  existed  near  the  toe  of  the  embankment 

and  could  spread  inward  into  the  center  of  the  embankment,  influ- 

* 

encing  overall  stability,  a  complete  nonlinear  analysis  of  the 
embankment  was  conducted. 

To  perform  the  nonlinear  analysis,  the  problem  was  con¬ 
ducted  in  two  phases,  the  first  being  concerned  with  the  dead 
load  stresses  developed  by  construction  of  the  embankment,  and 
the  second  phase  being  concerned  with  the  stresses  generated  by  • 
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the  storm  wave.  At  all  times  in  the  analyses,  total  intergranu¬ 
lar  stresses  were  maintained  to  assess  the  nonlinear  behavior  of 
the  soil  materials.  The  dead  load  solution,  in  addition,  was 
conducted  in  three  stages,  corresponding  to  the  potential  con¬ 
struction  loading  schedule  of  the  embankment.  A  portion  of  the 
finite  element  tuesh  used  for  these  analyses  is  shown  in  figure  58, 


6 . 1  Material  Propertie s 


For  the  various  soil  types  of  the  problem,  three  dif¬ 
ferent  soil  models  were  used  in  the  analyses.  For  the  embankment 
materials  (gravel,  stone,  dolose),  a  Coulomb-Mohr  material  law 
was  used  as  outlined  in  Section  3.3.  For  the  foundation  mate¬ 
rials  (sand  and  clay  strata)  nonlinear  soil  models  were  used 
which  produced  stiffening  under  increasing  confining  pressures 
and  softening  under  increasing  deviatoric  stresses  or  shear  load¬ 
ings.  These  soil  models  are  similar  to  the  nonlinear  models  of 
Section  3.4.  In  these  cases,  formulas  based  on  experimental  data 
(Ref.  13)  were  generated  for  typical  sand/clay  materials,  relat¬ 
ing  bulk  and  shear  moduli  for  the  soils  to  j'  ,  initial  void 

ratio,  etc.  For  the  clay  soil,  a  hyperbolic  stress-strain  model 
was  used. 


6*2  Dead  Load  Analysis 


The  dead  load  analysis  proceeds  in  the  usual  manner  as 
follows.  A  finite  element  mesh  is  superimposed  on  the  zone  of 
interests,  properties  of  each  element  specified  by  the  (in  gen¬ 
eral)  nonlinear  stress  strain  relations  of  the  materials  and 
suitable  boundary  conditions  applied  to  the  boundary  modes.  The 
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analysis  then  determines  the  solution  to  the  following  system  of 
equations. 


If*}  •  M{r\  -  {e*\ 


(47) 


where  h*\  is  the  vector  of  applied  loadings.  For  this  problem,  . 
the  applied  loading  is  simply  the  dead  loads  caused  by  the  mate¬ 
rial  weight.  The  vector  w  is  the  vector  of  node  displacements 
(unknown)  while  the  matrix  M  is  the  effective  stiffness  matrix 
of  the  finite  element  mesh  determined  from  the  initial  moduli  of 
the  materials.  The  vector  ^ F  |  are  correction  forceq  which  ac¬ 
count  for  the  fact  that  the  material  stress-strain  relations  vary 
with  applied  stresses  and  are  nonlinear.  To  solve  equations  (47) 
an  iterative  technique-  is  used  to  generate  the  nonlinear  correc- 

5 

tion  forces  with  node  displacement. 

The  zope  of  interest  consists  of  six  material  zones  and 
the  caisson.  The  caisson  was  simulated  by  including  the  finite 
element  mesh  through  the  structure  but  using  relatively  high 

V 

stiffness  properties  to  determine  rigid  body  response.  As  can  be 
noted  in  Fig.  57',  the  material  zone  was  considered  to  a  depth  of 
130  feet  below  M£iW,  through  the  second  relatively  stiff  founda¬ 
tion  sand  material.  For  the  static  analyses,  this  zone  cutoff 
was  considered  adequate  since  the  top  sand  and  day  layer  are 
relatively  soft  and  will  cause  the  majority  of  the  caisson  set¬ 
tlements.  In  addition,  from  the  previous  wave  pore  pressure  so¬ 
lution  (Section  5.0),  it  was  determined  that  wave  pore  pressure 
effects  will  not  extend  past  the  foundation  clay  layer.  Thus 
this  same  zone  definition  should  be  adequate  for  the  wave 
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liquefaction  analysis. 

\ 

The  lateral  extent  of  the  zone  of  interest  used  in  the 
analysis  is  larger  than  shown  in  Figure  58,  and  provided  a  foun¬ 
dation  extertt  of  about  200  feet  on  cither  side  of  the  embankment 
cross-section.  This  extent  was  deemed  adequate  to  allow  any  sta¬ 
bility  failures,  if  generated,  from  occurring.  The  boundary  con¬ 
ditions  on  the  lateral  boundaries  were  specified  as  vertical 
roller  supports  while  on  the  horizontal  boundary  they  were  speci¬ 
fied  as  hbrizontal  rollers. 

The  finite  element  mesh  consists  of  both  quadrilateral 
as  well  as  triangular  elements.  For  the  foundation  materials, 
single  element  thicknesses  were  provided  through  the  top  founda¬ 
tion  sand  and  clay  layers  while  two  elements  were  used  through 
the  bottom  foundation  sand  thickness.  For  this  study,  this  fine¬ 
ness  of  the  mesh  was  deemed  adequate  although  the  mesh  through 

the  clay  zone  should  probably  be  improved  in  more  refined 

* 

analyses. 


In  order  to  determine  settlement  characteristics  simu¬ 
lating  actual  constructive  procedures,  a  four-stage  analysis  was 
considered.  The  first  stage  was  concerned  with  determining  the 
initial  stress  and  strain  state  in  the  material  overburden  with¬ 
out  the  embankment.  This  was  accomplished  by  using  a  zero  unit 
weight  for  all  the  embankment  zones.  Once  this  solution  was 
available,  the  following  loadings  were  placed  on  the  structure  in 
the  sequence  that  could  be  used  during  actual  construction.  Dur¬ 
ing  each  of  the  four-stage  analyses,  a  double  iteration  procedure 
was  used  to  improve  the  accuracy  of  the  results.  Applied  load¬ 
ings  during  each  stage  were  applied  in  three  increments  and 
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numerical  interation  (equation  47)  was  used  during  each  load 
increment. 

The  first  construction  stage  loading  considered  the  ap¬ 
plication  of  the  loads  caused  by  the  caisson  plus  a  15-foot  thick 
layer  of  the  fill  material  as  shown  in  Figure  59.  The  factors  of 
safety  computed  in  each  element  are  also  shown  in  the  figure. 

The  factor  of  safety  is  defined  as  shear  strength  available  (due 
to  the  hydrostatic  compression  stress)  divided  by  the  effective 
shear  stress  applied.  As  may  be  noted  from  Figure  59,  the  fac¬ 
tors  of  safety  throughout  the  zone  are  adequate  (greater  than 
1.0)  except  in  the  gravel  material  immediately  below  the  caison. 
This  indicates  that  the  properties  assumed  for  this  material  are 

t 

too  soft  to  adequately  restrain  the  caisson.  A  factor  of  safety 
of  1.0  indicates  that  the  stress  state  is  on  the  yield  surface 
and  relatively  large  deformations  are  possible.  Preliminary  in¬ 
vestigation  indicated  that  the  properties  used  in  the  analyses 
are  conservative  and  higher  stiffness  and  strength  properties 
could  be  used.  However,  this  also  implies  that  the  placement  of 
this  material  must  be  such  as  to  ensure  the  adequacy  of  the 
properties. 

The  second  construction  stage  loading  is  shown  in 
Figure  60  and  the  results  are  similar  to  the  above,  namely  that 
the  safety  factors  in  the  gravel  are  about  1.0.  The  final  con¬ 
struction  stage  loading  is  shown  in  Figure  61  together  with  all 
the  safety  factors  throughout  the  embankment.  As  may  be  noted, 
the  conservative  property  estimates  for  the  gravel  fill  material 
lead  to  low  factors  of  safety  throughout  the  zone. 


The  settlement  profile  of  the  caisson  during  each 
loading  stage  is  shown  in  Figure  62.  The  results  indicated  that 
the  elastic"  properties  assumed  for  the  caisson  were  adequate 
in  that  the  caisson  moved  essentially  as  a  rigid  body.  As  can 
be  seen,  the  settlements  predicted  for  the  caisson  vary  from 
7.5  inches  on  the  inside  (left  side  of  the  mesh)  to  14.6  inches 
on  the  side  toward  the  embankment  centerline.  This  large  settle¬ 
ment  is  due  primarily  to  the  settlements  caused  by  the  shear  flow 
indicated  by  the  low  safety  factors  in  the  gravel  materials. 

Again,  using  stiffer  values  for  this  data  will  lead  to  correspond¬ 
ingly  smaller  caisson  motions  but  placement  conditions  must  be 
such  so  as  to  guarantee  the  availability  of  the  properties. 

•  i 

6*3  Wave  Liquefaction  Analysis 

Prior  to  starting  the  liquefaction  analysis,  several 

modifications  were  mace  to  th  i  intergranular  soil  properties  of 

the  embankment  fill  materials.  The  initial  elastic  moduli  of 

these  materials  were  increased  by  50%  while  the  friction  angl'ri 

for  the  quarry  run  material  was  increased  from  40°  to  45°.  This 

was  done  since  the  safety  factors  in  the  gravel  material  were 

1.0  (stress  points  on  the  yield  surface)  and  initial  liquefaction 

analyses  indicated  that  the  embankment  could  not  withstand  any 

!  » 

& ^-9-^ f icsnt  wsvg  motion*  To  propGjrly  mclud6  this  GffGct,  thG 
dead  load  solution  should  be  rerun,  to  obtain  a  consistent  set 
of  initial  conditions  for  the  liquefaction  analysis.  However, 
for  this  demonstration  investigation,  it  was  felt  that  this 
improvement  was  unnecessary. 

For  this  embankment  problem,  pore  pressure  time 
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histories  were  specified  along  the  top  of  the  foundation  sand 
since  it  was  felt  that  the  permeability  effects  through  the 
coarse  embankment  materials  were  sufficiently  nonlinear  as  to 
preclude  a  realistic  analysis  based  upon  D'arcy's  law.  Thus 
seepage  flow  assumed  to  occur  within  the  foundation  soils  only 
and  not  through  the  embankment  materials.  Excess  pore-pressure 
histories  were  specified  at  the  locations  shown  in  Figure  57. 
These  pressures  were  specified  by  essentially  taking  the  storm 
wave  pressures  and  scaling  them  linearly  with  distance  into  the 
embankment  and  adding  a  small  time  lag  to  phase  the  pressures. 

It  is  clear  that  to  properly  analyze  this  problem,  experimental 
data  would  be  required  to  ascertain  a  complete  set  of  pore 
pressure  boundary  conditions. 

The  problem  was  started  with  the  wave  applied  upstream 
of  the  embankment  (right  side  of  the  mesh)  and  moved  into  the 
embankment  in  relatively  small  time  steps.  The  numerical  results 
shown  are  referenced  to  time  zero  occurring  when  the  wave  crest 
is  approximately  above  the  toe  of  the  embankment.  Figures  63 
through  70  indicate  the  factors  of  safety  computed  in  each  ele¬ 
ment  as  the  wave  progresses  into  the  embankment.  Numerical  sta¬ 
bility  problems  appear  to  begin  when  the  negative  pore  pressures 
are  generated  on  the  right  side  of  the  mesh  (above  the  toe) . 

The  displacements  computed  indicate  the  clay  layer  under  the 
embankment  moving  down  and  to  the  right  with  the  material  to 
the  right  of  the  toe  moving  upward,  causing  a  general  decrease 
in  the  factors  of  safety  in  the  embankment  materials.  Further 
computation  tends  to  exaggerate  these  results  leading  to  an 
eventual  collapse  of  the  embankment.  Thus  the  calculations 


indicate  an  erosion  of  the  foundation  material  strengths  under  the 
toe  of  the  embankment.  This  is  caused  by  the  seepage  pressures 
generated  from  the  water  flow  moving  from  under  the  embankment 
(positive  pore  pressures)  . 
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7.0  SUMMARY 


This  report  has  presented  the  results  of  a  study  to  in¬ 
clude  the  effects  of  pore  water  on  the  response  of  soils  to  tran¬ 
sient  loadings.  Numerical  results  have  been  presented  for  prob¬ 
lems  which  fall  into  different  categories,  in  an  attempt  to  span 
the  problem  areas  of  general  interest.  The  first  set  of  problems 
was  concerned  with  computation  of  soil  response  to  loadings  in 
configurations  which  would  be  encountered  in  the  laboratory, 
namely,  triaxial  and  uniaxial  configurations. 

The  second  set  of  problems  was  concerned  with  the  im¬ 
pact  of  entrapped  pore  fluid  on  wave  propagation  velocities.  As 
goes  without  saying,  the  wave  motion  is  completely  governed  by 

i 

the  boundary  conditions  since  the  fluid  has  been  assumed  incom¬ 
pressible.  For  the  problem  conditions  assumed,  the  wave  speeds 
in  the  soil  are  much  smaller  than  the  actual  wave  speeds  in  water 
so  that  the  incompressible  assumption  for  the  water  appears  rea¬ 
sonable.  For  the  case  where  soil  wave  speeds  approach  that  of 
the  entrapped  fluid,  this  assumption  becomes  questionable  and  at 
least  fluid  compressibility  must  be  included  in  the  analysis. 

The  last  two  problem  sets  have  to  do  with  a  study  of 

submerged,  off-shore  soils  to  transient  ocean  storm  waves,  a 

» 

problem  area  receiving  more  and  more  consideration  as  off-shore 
construction  increases  (platforms,  nuclear  reactor  stations,  pro¬ 
tective  dikes,  etc.).  The  approaches  used  appear  to  offer  an  at¬ 
tractive  solution  technique  to  this  important  class  of  problems. 

One  aspect  of  the  numerical  computation  problem  which 
has  not  been  discussed  in  this  report  concerns  the  stability  of 


the  computations.  It  is  quite  clear  that  such  .considerations 
were  made  during  the  course  of  the  study  to  arrive  at  an  accept¬ 
able  computation,  but  no  formal  study  of  the  problem  was  conducted 
during  this  project.  The  stability  considerations  have  to  do 
with  the  types  of  pioblems  which  may  be  studied  as  well  as  the 
time  step  required  for  the  transient  analyses  to  achieve  a  suc¬ 
cessful  integration  procedure.  Some  information  has  appeared  in 
the  literature  concerning  the  first  problem  (Refs.  9  and  10) 
while  no  information  on  the  second  has  been  seen.  Thus  approxi¬ 
mate  rules  of  thumb  have  been  devised  for  specific  problem  types 
but  are  not  general  enough  as  yet  to  be  considered  adequate. 
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A.  1 


APPENDIX  A 


FORMULATION  OF  SYSTEM  EQUATIONS 


In  the  following  presentation,  the  analysis  will  be  car¬ 
ried  forth  for  a  typical  element  of  the  free-field  mesh.  The  dis¬ 
placement  field  for  the  element  is  assumed  to  be  linear  and  the 
displacement  of  any  point  within  the  element  can  be  written  as 


v  *  bl*  W 


( A .  1 ) 


where  (U,w*)  are  the  horizontal  and  vertical  displacement  components 
and  ^  are  each  a  set  of  arbitrary  coefficients,  with  the 

number  of  coefficients  equal  to  the  number  of  element  vertices  to 
provide  the  proper  number  of  degrees  of  freedom  for  the  element. 

The  vector  is  formed  by  a  proper  set  of  element  func¬ 

tions  and  depend  upon  the  element  type  being  considered.  For  a 
typical  triangular  element  (Fig.  A.l)  this  vector  is 


bV  b-r>b 


(A. 2) 


while  for  a  typical  rectangular  element  (Fig.  A. 2) 


(A.  3‘) 
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A. 2 

For  a  general  quadrilateral  element  (Fig.  A.  3)  the  vector  spec 
ified  by  equation  A. 3  is  used  m  the  transformed  coordinate  system 
By  substituting  the  coordinates  of  the  nodes  into  equation  A.l, 
the  coefficients  ^  ^  can  be  replaced  as  unknowns  by  the  node 

point  displacement  components ,  or 


(A. 4) 


This  simple  displacement  function  assumed  for  the  element  allows 

i 

for  determining  any  interior  displacement  in  terms  of  the  nodal 
displacements  and  ensures  that  the  displacements  between  any  two 
adjacent  elements  will  be  continuous  for  any  arbitrary  specifica¬ 
tion  of  nodal  displacements.  Higher  order  element  formulations 
are  also  available  to  satisfy  the  above  criteria. 

The  strains  developed  at  any  point  within  the  element 
can  be  determined  from  the  strain  displacement  relations  for  the 
particular  configuration,  or 

(A*5) 

where  is  the  strain  vector  with  components 

i€r,€J\  6) 


The  superscript  T  in  equation  A. 6  indicates  total  strains 


A.  3 


For  the  combined  stress/pore  pr.  ss  pro:  lero,  r to  intergranular 
stress  are  related  to  the  por  ■  pro  ,.uro  by 


(crj  *  {,f\ -  nr \  £ \  (A.7) 

where  ^(r|  are  the  total  stress  in  the  body  defined  by 

i0"]  *{o*r»  ,  Cj ,  rrj|  (A. 8) 

w)  are  the  effective  or  intergranular  stresses  and  T  is  the  pore 

•  I 

pressure.  The  vector  is  defined  as  Jl ,  1  ,  1 ,  Of.  The  effec¬ 
tive  stresses  are  related  to  the  strains  through  the  general  stress 

strain  relations 


W*Lc]l«T-«") 


where  are  defined  as  the  nonlinear  components  of  the  total 

s^rain  an<*  [C]  is  the  usual  elastic  sties  .  strain  matrix  which,  for 
example,  can  be  defined  for  the  axisyrtunetnc  problem  by 
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(A. 10) 


and 


where  E  is  Young's  modulus  and  2>  is  Poisson's  ratio.  For  relatively 
simple  material  models  (such  as  The  Mises  or  Coulomb-Mohr  plastic 
models),  the  nonlinear  strains  represent  the  nonrecoverabi "  or 
plastic  strain  components.  For  more  complicated  material  models, 
the  vector  represents  a  fictitious  set  of  strains  required  to 

yield  the  proper  stresses. 

To  satisfy  equilibrium  conditions  at  the  element  nodes 
with  the  total  stress  field  within  the  element,  the  usual  virtual 
work  principal  is  used.  The  internal  work  performed  by  the  stresses 
on  a  virtual  total  strain  field  is  defined  by 

JW;  *.  J  (£<=T}Vl  M  (a. in 

where  the  integral  is  taken  over  the  element  volume.  The  corre¬ 
sponding  external  work  performed  by  forces  applied  at  the  nodes  is 

£*We  s  (A. 12) 

where  are  the  horizontal  force  components  at  each  node  and 

are  the  corresponding  vertical  force  components.-  Equating 

i 

the  internal  and  external  work  expressions  and  making  us  of  the 
definitions  previously  described,  the  force  components  that  must 
be  applied  at  each  node  point  to  maintain  equilibrium  w ith  the 
total  stress  within  an  element  are 


fa  Jr  t  ♦  faoJlfa*  -faJ\  -{*V\ 
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(A. 13) 


The  matrices  t  Lte.tft)^  >  etc.  are  the  usual  elastic  stiffness 

matrices  and  are  defined  by 


t  l!jV ■  UBiVtcUBj)  W 


(A. 14) 


where  the  subscripts  (i,j)  take  on  the  values  of  (u,w)  .  The  terms 
fa '  faurl  represent  the  correction  forces  to  account  for  mate¬ 
rial  nonlinearity  and  are  defined  by 


(A. 15) 


where  again  the  subscript  (i)  takes  on  the  values  (u,w) .  The  terms 
fa*} '  represent  the  effects  of  pore  pressure  on  the  equilib¬ 
rium  equations  and  are  defined  by 


’  I  TT  AV 


(A. 16) 


In  the  computer  program  developed,  the  pore  pressure  variation  is 
assumed  to  be  a  linear  one  over  the  element  or 


T-lcjV  hVt-tQU] 


(A. 17) 


A. 6 


where  the  vector  |TiTj  represents  th>  nodal  point  pore  pressures. 
Substituting  equation  A. 17  into  A.!6  then  yields 


un  -Lkiiw 


(A. 18) 


where  (. <W 


(A. 19) 


where  the  subscript  i  represents  both  the  u  and  w  directions. 

To  relate  the  pore  presusres  to  the  node  point  displace¬ 
ments  (or  velocities),  seepage  effects  are  considered.  The  seepage 

\ 

equations  are  obtained  by  minimizing  the  functional  (Ref.  8). 


11  -  L  U  [  ^  ♦  2  -  TT  QW 


(A. 20) 


The  permeability  components  krr ,  kr^  »  k^  are  related  to  the  prin¬ 
cipal  permeability  coefficients  by 

—  k|  5^  Q  i  is  j_s.cS  £> 

kr*  *  SvAfttfes.©  Cfei-k.*)  (A>21) 

i 

k^  *  k.cos'©  + 


where  (k|f  kt.)  are  the  principal  coefficients  in  two  orthogonal 
directions  and  S  is  the  angle  from  the  i  -direction  to  the  f -direction 
positive  in  the  clockwise  sense. 

For  a  particular  element,  the  variation  of  pore  pressure 
over  the  element  is  assumed  to  be  linear  or 


A. 7 


(A. 22) 


Substituting  equation  A. 22  inva  A. 20,  the  functional  can  be  written 
as 


A  *i  M'M WLo]  \t} .  (a. 23) 

where  the  matrix  ^0  is  defined  by 


The  vector  is  defined  as  the  volume  integrals 

1T1  *  Iv  M  (A. 25) 

where  Q  is  the  volume  decrease  per  unit  volume  per  unit  time. 

The  solution  to  the  seepage  problem  is  obtained  by  mini¬ 
mizing  the  function il  with  respect  to  the  nodal  pressures  leading 
to  the  conservation  equation 

r 

faVlS]fa]  «■  (a.  26) 

or  to  further  condense  notation 


LhH.it}  ■  [oVL-j'l 


v A . 27 ) 


A. 6 


The  vector  can  be  evaluated  by  considering  the  volume 
compression  of  the  solids  plus  that  of  the  water,  or 


*  -  -  WUl 


(A. 28) 


where  E  v  is  the  effective  bulk,  modulus  of  water  and  the  vector 
r  *  t) 

Jis  the  total  strain  rate  vector.  The  strain  rates  are  related 
to  the  nodal  velocities  by 


[*T]  8  LfculW  *  iBw/ftur} 

% 


Substituting  the  above  into  the  conservation  equation  leads  to 


WM LfejV\  -  w-}  +  jyiw 

[H¥h  Lt>TLt3\_t>1 


(A. 30) 


(a)  Original  Coordinates 


Co)  Transformed  CoerJinates 

Fig.  A.  3  .General  Quadrilateral  Elc^'  ’  " 
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